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This  final  technical  report  documents  the  research 
efforts  in  the  numerical  simulations  of  tropical  weather 
systems  carried  out  under  the  auspices  of  ONR  Contract 
N0001 4-82-C-2306,  for  the  Atmospheric  Physics  Branch, 
Space  Science  Division,  of  the  Naval  Research  Laboratory 
(NRL) . 


Our  research  efforts  mainly  concentrate  in  two  major 
study  areas:  the  interactions  between  two  tropical 
cyclones  and  the  construction  of  an  axisymmetric  ocean 
model.  The  research  in  both  of  these  areas  has  been 
completed  and  reported  in  two  articles.  A  paper  entitled 
"A  Numerical  Study  of  the  Interactions  Between  Two 
Tropical  Cyclones"  was  accepted  by  the  Monthly  Weather 
Review  and  Is  to  be  published  in  the  1983  September  issue. 
The  draft  of  a  second  article  on  the  ocean  model  has  been 
completed.  Both  articles  are  included  as  Appendices  in 
this  final  report.  We  have  also  attached  a  listing  of  the 
computer  code  of  the  ocean  model,  which  is  also  stored  in 
and  accessible  from  the  TI-ASC  of  NRL.  Here,  a  brief 
summary  of  the  study  result  will  be  given. 

The  interactions  between  atmospheric  vortex  pairs  are 
simulated  and  studied  with  a  nondivergent  barotropic  model 
and  a  three-dimensional  tropical  cyclone  model  (NRL/SAI 
mesoscale  model).  Numerical  experiments  with  nondivergent 
barotropic  vortex  pairs  show  that  the  relative  movements 
of  the  vortices  are  sensitive  to  the  separation  distance 
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and  the  characteristics  of  the  swirling  wind  of  the 
vortex.  No  observed  mutual  attraction  is  found  in  any  of 
the  nondivergent,  barotropic  vortex  pairs  tested. 

Results  from  the  3D  NRL/SAI  tropical  cyclone  model 
show  that  on  a  constant-f  plane  with  no  mean  wind,  the 
movements  of  the  two  interacting  tropical  cyclones  consist 
of  a  mutual  cyclonic  rotation,  attraction,  and  eventual 
merging,  in  agreement  with  Fujiwhara's  description.  The 
displacement  of  one  interacting  storm  in  the  mutual 
rotation  is  proportional  to  the  combined  strength  of  the 
binary  system,  but  inversly  proportional  to  the  size  of 
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the  storm  and  to  the  square  of  the  separation  distance. 
The  rate  of  merging  is  related  to  the  development  of  a 
mean  secondary  circulation  on  the  radial-vertical  plane, 
and  is  quite  independent  of  the  strength  of  the  two 
tropical  cyclones. 

The  latitudinal  variation  of  the  Coroilis  parameter 
adds  a  northwest  beta  drift  to  the  trajectories. 
Depending  on  their  relative  strength  and  location,  the 
beta  drift  can  either  speed  up  the  merging  process  or 
separate  the  two  interacting  tropical  cyclones. 

The  axisymmetric  ocean  model  consists  of  primitive 
equations  for  the  conservation  of  momenta  in  three  spatial 
dimensions  and  the  buoyancy.  A  Boussineq  assumption  is 
made  so  that  the  background  stratification  is  kept 
constant,  the  horizontal  and  vertical  diffusion  is  of  the 
Fickian  type. 


A  leapfrog  temporal  integration  is  employed.  The 
grid  is  fully  staggered  as  Arakawa  C  type.  The  system  is 
non-hydrostatic,  the  resultant  elliptic  equation  for  the 
pressure  is  solved  by  a  stablized  error  vector  propagation 
technique.  The  basic  equations,  the  finite  differencing 
form,  and  boundary  conditions  are  discussed  in  detail  in 
the  attached  Appendix. 


APPENDIX  I 

A  NUMERICAL  STUDY  OF  THE  INTERACTIONS  BETWEEN 
TWO  TROPICAL  CYCLONES 
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ABSTRACT 


The  interactions  between  atmospheric  vortex  pairs  are  simulated  and 
studied  with  a  nondivergent  barotropic  model  and  a  three-dimensional  tropical 
cyclone  model. 

Numerical  experiments  with  nondivergent  barotropic  vortex  pairs  show 
that  the  relative  movements  of  the  vortices  are  sensitive  to  the  separation 
distance  and  the  characteristics  of  the  swirling  wind  of  the  vortex.  No 
mutual  attraction  is  found  in  any  of  the  nondivergent,  barotropic  vortex  pairs 
tested. 

Results <from  the  3D  tropical  cyclone  model  show  that  on  a  constant- f 
plane  with  no  mean  wind,  the  movements  of  the  two  interacting  tropical  cyclones 
consist  of  a  mutual  cyclonic  rotation,  attraction,  and  eventual  merging,  in 
agreement  with  Fujiwhara'a  description.  The  displacement  of  one  interacting 
storm  in  the  mutual  rotation  is  proportional  to  the  combined  strength  of  the 
binary  system,  but  inversly  proportional  to  the  size  of  the  storm  and  to  the 
square  of  the  separation  distance.  The  rate  of  merging  is  related  to  the 
development  of  a  mean  secondary  circulation  on  the  radial-vertical  plane,  and 
is  quite  independent  of  the  strength  of  the  two  tropical  cyclones. 

The  latitudinal  variation  of  the  Coriolis  parameter  adds  a  northwest 
beta  drift  to  the  trajectories.  Depending  on  their  relative  strength  and 
location,  the  beta  drift  either  speeds  up  the  merging  process  or  separates 
the  two  Interacting  tropical  cyclones^ 
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1.  INTRODUCTION 

When  two  tropical  cyclones  are  present  simultaneously  in  the  same 

region,  it  is  often  observed  that  they  rotate  around  each  other  with  decreasing 

separation  between  them  in  the  absence  of  large  scale  wind  flow  (Fig.  1).  The 
phenomenon  was  made  well-known  by  Fujiwhara  (1921),  and  is  therefore  referred  to 

as  the  Fujiwhara  effect.  By  laboratory  experiment  and  geophysical  observation, 
Fujiwhara  (1923,  1931)  demonstrated  that  the  relative  motion  of  two  counterclock¬ 
wise  vortices  was  a  counterclockwise  rotation.  Haurwitz  (1951)  examined  several 
tropical  cyclone  pairs  by  introducing  the  concept  of  center  of  mass  around  which 
the  two  tropical  cyclones  rotate  about  each  other.  By  approximating  the  circu¬ 
lation  around  a  tropical  cyclone  with  that  of  a  Rankine  vortex,  Haurwitz  (1951) 
derived  a  relationship  between  the  rotation  rate  and  the  sum  of  the  total  mass 
circulation  of  the  two  tropical  cyclones.  Many  discrepancies  were  found  when 
he  applied  the  relationship  to  observations.  Haurwitz  attributed  the  discrepan¬ 
cies  to  the  Influence  of  large  scale  flow  and  lack  of  data,  which  led  to  deficiency 
in  analyses. 

Hoover  (1961)  studied  binary  tropical  cyclones  in  both  the  Atlantic  and 
Western  Pacific  Oceans.  He  found  that  the  interaction  between  tropical  cyclone 
pairs  in  the  Western  Pacific  Ocean  agrees  with  Fujiwhara *s  description  while 
those  pairs  in  the  Atlantic  Ocean  rotated  in  an  anticydonic  sense.  He 
suggested  that  the  different  large  scale  atmospheric  flow  patterns  in  the 
two  basins  may  have  caused  the  binary  systems  to  behave  differently.  The  in¬ 
fluence  of  the  large  scale  flow  was  also  noted  by  Liu  and  Wang  (1966).  They 
found  that  two  interacting  tropical  cyclones  in  the  Western  Pacific  are  not 
always  attracted  to  each  other  when  there  are  strong  shears  in  the  environmental 
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flow.  Recently,  Dong  and  Neumann  (1982)*  found  that  storm  pairs  exhibiting 
behavior  most  in  accordance  with  Fujiwhara's  description  were  located  in 
the  Intertropical  Convergence  Zone  where  horizontal  shears  in  large  scale 
flow  are  negligible.  They  suggest  that  the  effects  of  environmental  flows 
be  filtered  before  the  real  Fujiwhara  effects  can  be  determined.  But  to  define 
and  remove  the  large  scale  flows  from  observational  data  is  difficult  to 
accomplish. 

Over  the  35  year  period  1946-1981,  storm  pairs  known  to  have  interactions 
averaged  1.5  annually  in  Western  North  Pacific  and  0.33  annually  in  the  Atlantic 
(Dong  and  Neumann,  1982).  The  presence  of  binary  interacting  has  been  noted  to 
have  contributed  to  forecast  errors  of  tropical  cyclone  tracks  (Brand,  1970; 
Jarrell  et.  al. ;  1978;  Neumann,  1982).  Forecasting  as  well  as  analyzing  a 
single  tropical  cyclone  is  often  hindered  by  the  paucity  of  observational  data 
in  the  tropical  cyclone  basin  (Neumann,  1982);  the  presence  of  two  storms  in 
close  proximity  can  further  compound  the  difficulties. 

The  purpose  of  our  study  is  to  investigate  the  interactions  between  two 
tropical  cyclones  by  numerical  simulations.  Because  the  spatial  resolutions  of 
the  models  are  better  than  the  current  observational  network,  and  because  nu¬ 
merical  models  can  be  controlled  to  produce  "clean"  results  void  of  undesirable 
factors,  analyzing  realistic  numerical  simulations  can  sometimes  result  in  a 
better  isolation  and  understanding  of  the  phenomenon  than  can  be  acheived  from 
an  observational  approach.  In  this  paper  we  will  first  determine  the  role  of 
vorticity  advection  between  the  two  vortices.  For  this  purpose,  a  nondivergent, 
barotropic  model  is  introduced  to  test  two  types  of  vortex  pairs  with  different 

*Dong,  K.  and  C.  J.  Neumann,  1982:  On  the  relative  motion  of  binary  tropical 
cyclones.  Regional  Scientific  Conference  on  Tropical  Meteorology,  Tsukuba 
Ibaraki,  Japan,  Oct.  1982. 

•  2 


swirling,  winds .  These  barotropic  tests  will  be  presented  in  Section  2. 

In  Section  3,  three-dimensional  simulations  of  the  interactions  between  two 
diabatically  driven  tropical  cyclones  on  a  constant-f  plane  and  with  variable 
f  will  be  discussed.  Our  findings  will  be  summarized  in  Section  4. 


2.  INTERACTIONS  BETWEEN  NONDIVERGENT ,  BAROTROPIC  VORTEX  PAIRS 

In  this  section,  we  investigate  the  interactions  between  nondivergent 
barotropic  vortex  pairs.  Through  the  interactions  of  such  vortex  pairs,  we  can 
determine  the  contribution  of  horizontal  advection  of  vorticity,  because  in 
such  a  system  advection  is  the  only  mechanism  for  interaction.  A  description 
of  the  nondivergent,  barotropic  model  will  be  presented  first,  and  the  experi¬ 
mental  design  and  the  results  will  then  be  discussed. 


a.  Nondivergent  Barotropic  Model 


The  simple  non-divergent,  barotropic  model  can  be  described  as 

3  2  2 

-  V  i|>  -  -V^,  .  (V  <{>  +  f) ,  and  i 
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-  fc  x  2^  (2) 

where  f  is  the  Coriolis  parameter,  is  the  stream  function,  Vib  is  the 

nondivergent  wind,  and  ft  is  a  vertically  pointing  unit  vector.  The  bound 


where  the  relative  vorticity  is  defined  by 


(4) 
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is  solved  by  a  stabilized  error  vector  propagation  method  (Iladala,  1978). 


b.  Experimental  Design 

The  major  application  of  the  nondivergent  barotropic  model  is  to  determine 
the  effects  of  separation  distance  and  the  radial  distribution  of  tangential 
winds  on  the  interaction  of  the  two  vortices.  Two  kinds  of  wind  distributions 
were  tested.  The  first  kind  (type  A)  of  vortex  is  defined  by  its  cyclonic 
swirl  wind  vQ  as  function  of  radius  r  from  the  vortex  center 

tfr~"  '  •  -  ..  __ 

Ar  (1  -  sin  - ) ,  0  <  r  <  r«  ' 

r  -  -  o  ,  . 

o 


vo  - 


(5) 


0,  otherwise 

where  constant  A  -  4  x  IQ”4  s-i  and  rQ  =  400  km.  Equation  (5)  yields  a 
Tna-r-rtnirm  swirl  of  ^26  m  s-*  at  r  ■  150  km  and  a  maximum  vorticity 
5  I  7,2  x  10“4  s*^  at  r  »  0.  We  note  that  there  is  a  cutoff  of  vo  at  rQ. 
The  second  kind  of  vortex  (type  B)  is  defined  as 

r  ^ 

'  v  -  Br  exp  ( - 

°  *; 


(6) 


where  the  e-folding  distance  rg  is  150  km.  By  letting  the  constant 
B  _  3.6  x  10-*  s'1,  (6)  yields  a  vortex  with  similar  strength  as  that  described  by 
(5)  with  maximum  swirl  of  ^29  m  s'1  and  a  maximum  C  =  7  x  10"4  s-1.  Type  B  vortex 
differs  from  Type  A  in  that  there  is  no  cutoff  of  swirl.  Fig.  2  compares  the  radial 
distributions  of  relative  vorticities  described  by  (5)  and  (6). 
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Four  initial  separation  distances  (300,  400,  600,  and  1000  km)  have 
been  tested  for  each  type  of  vortex  pairs.  All  intergrations  with  the  baro- 
tropic  model  are  performed  with  constant  f  ■  4.37  x  10”^  s~*. 

c.  Results 

Fig.  3  shows  the  trajectories  of  storm  pairs  having  two  types  of  swirl 
wind  at  four  separation  distances.  It  is  very  clear  from  Fig.  3  that  the 
smaller  the  separation  distance  the  faster  the  mutual  transport.  For  in¬ 
stance,  at  a  separation  distance  of  1000  km,  neither  type  A  vortex  pair  (with 
swirl  cutoff  at  r  *  300  km)  nor  type  B  vortex  pair  can  induce  mutual  motion. 

But  at  a  separation  distance  of  400  km,  they  move  at  a  speed  of  ^<400  km  day-* . 

It  is  also  evident  that  the  mutually-induced  motions  of  type  A  and  type  B 
vortices  are  very  different,  in  spite  of  the  values  of  constants  for  A  and  B 
which  were  chosen  to  give  vortices  of  similar  strength.  Furthermore,  the 
trajectories  of  type  A  vortex  pairs  are  more  anticydonic.  This  may  be  a 
result  of  the  fact  that  type  B  vortices  have  positive  vorticities  at  r  <  200  km 
whereas  the  vorticities  of  type  A  vortices  change  sign  at  r  *  150  km  (Fig.  2). 
Only  vortex  pairs  at  small  separation  distances  rotate  in  a  cyclonic  fashion 
because  they  interact  with  positive  shears.  The  motion  of  vortices  in  our 
model  can  only  be  caused  by  the  advection  of  vorticity,  the  shear  in  one  vortex 
can  very  much  determine  the  movement  of  the  other.  These  results  indicate  that 
the  mutual  motion  of  two  interacting  nondivergent ,  barotropic  vortex  pair  are 
quite  sensitive  to  the  charaterlstics  of  the  swirl  winds. 

In  all  the  experiments  illustrated  in  Fig.  3  the  storm  pairs  drift  apart, 
there  is  no  mutual  attraction  as  observed  in  some  interacting  typhoons.  This 
suggests  that  the  observed  mutual  attraction  in  typhoon  pairs  may  be  due  to 


the  divergence  and/or  convergence  that  is  not  included  in  the  barotropic 
model.  Indeed,  complicated  diabatic  processes  in  tropical  cyclone  such  as 
long  wave  radiation,  surface  boundary  layer  effects  and  moist  convection 
generate  convergent  flow  in  lower  troposphere  and  divergent  flow  in  upper 
troposphere.  The  irrotational  component  of  the  vortex  circulation  may  be 
responsible  for  the  occurrence  of  the  observed  cyclonic  rotation  and  mutual 
attraction. 


3.  INTERACTIONS  BETWEEN  TROPICAL  CYCLONE  PAIRS 

We  have  seen  that  nondivergent  pairs  do  not  cause  a  mutual  motion 
similar  to  the  description  of  Fujiwhara.  The  observed  Fujiwhara  effects 
may  be  due  to  dynamics  that  can  only  be  resolved  by  a  more  complete  model. 
To  see  that,  we  will  simulate  the  interactions  between  two  diabatically 
driven  tropical  cyclones  with  a  three-dimensional  model. 


a.  Three-dimensional  Tropical  Cyclone  Model 


The  baroclinlc  model  is  identical  to  the  one  in  Chang  and  Madala  (1980) 
and  Chang  (1982),  except  for  parameterization  of  the  latent  heating.  The 
governing  equations  are  in  surface-pressure-weighted  flux  form  for  conser¬ 
vation  of  momentum,  temperature  and  water  vapor.  The  normalized  pressure 
o  ■  p/p„  is  the  vertical  coordinate,  where  p  is  the  surface  pressure.  The 
system  is  assumed  hydrostatic.  The  bulk  boundary  layer  parameterization  is 
based  on  a  generalized  similarity  theory  (Chang,  1981).  The  model  has 
51  x  51  horizontal  grid  points  with  seven  sigma  layers  in  the  vertical.  The 
horizontal  resolution  if  %°  in  both  the  latitudinal  and  longtitudinal  directions 
The  east-west  boundaries  are  cyclic.  The  boundary  conditions  at  the  north  and 
south  boundaries  are  such  that  the  second  derivatives  of  thermodynamic  variables 
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normal  to  the  boundaries  vanish.  In  addition,  diffusion  coefficients  are 
increased  near  the  north  and  south  boundaries  to  damp  numerical  noise  there. 

Kuo's  parameterization  was  used  in  Chang  and  Madala  (1980)  and  Chang 
(1982),  but  a  prescribed  heating  is  applied  in  this  model  as  done  by  Anthes 
(1971)  In  a  axisymmetric  model.  The  heating  rate  here  is  defined  as 


Q(r,a) 


nr  —  . 

4o  eos(— )  .in  0(0-0,  i)J,  f„  r<R>  o.l  jo  <0.1 


0,  otherwise 


where  r  is  the  distance  between  a  grid  point  and  the  low  pressure  center, 
and  R  -  300  km  is  the  limit  of  the  heating  function.  Two  values  of  $Q,  100  K  day”" 
and  200  K  day-*,  have  been  used  in  various  numerical  experiments  to  define  the 
weak  and  strong  tropical  cyclones,  respectively.  The  vertical  and  horizontal 
distributions  of  the  heating  pattern  described  in  (7)  are  illustrated  in  Figs. 4 
and  5.  The  vertical  heating  distribution  is  similar  to  that  of  the  differences 
between  temperatures  in  convective  clouds  (Tc)  and  the  environment  (T)  in  a 
mean  hurricane  season  sounding  for  the  Gulf  of  Mexico  as  computed  by  a  one¬ 
dimensional  cloud  model  (Anthes,  1977,  Fig  4a).  The  horizontal  heating  dis¬ 
tribution  agrees  with  the  mean  rainfall  rate  Inferred  from  satellite  observa¬ 
tion  in  a  typhoon  (Adler  and  Rogers,  1977),  except  for  the  observed  smooth 
fall-off  at  r  >  300  km.  No  effort  is  attempted  to  simulate  the  eye  because 
of  the  model  horizontal  resolution. 

The  heating  prescribed  by  (7)  nevertheless  generates  realistic  circu¬ 
lations  for  tropical  cyclones.  Figure  6  shows  the  radial  distribution  of 
the  quasi-steady  wind  speeds  at  the  sixth  (a  *  0.85)  and  seventh  (o  •  0.965) 
model  layers  after  24  h  of  heating  with  <$0  *  200  K  day-*.  The  wind  speeds 
have  a  peak  at  r  ■  1°  and  decrease  gradually  outward  without  discontinuity 


at  r  ■  R  *  300  km. 

We  note  however  that  by  using  the  prescribed  heating  in  (7)  the  effects 
of  the  interaction  between  the  two  cyclones  on  the  scale  of  the  cumulus  con¬ 
vection  cannot  be  adequately  simulated.  In  reality  the  momentum  field  in 
each  storm,  which  affects  the  cumulus  convection,  can  be  modified  by  the 
proximity  of  another  storm.  A  change  in  the  cumulus  convection  in  each  storm 
may  alter  the  storm's  intensities,  which  can  in  turn  affect  the  interaction 
between  the  two  tropical  cyclones.  But  these  feedbacks  may  be  secondary  and 
are  only  important  when  the  separation  of  the  two  storms  is  small.  As  a  pre¬ 
liminary  study,  the  more  economical,  prescribed  heating  is  used  to  investigate 
the  first  order  effects  in  the  interaction. 


b.  Experimental  Design 

The  tropical  cyclone  pairs  in  all  numerical  experiments  (Table  1)  with 
the  3D  model  are  dynamically  Initialized  by  a  24h  stationary  heating  at  two 
locations,,  i.e.,  by  applying  (7)  at  two  fixed  grid  points  for  24h.  Dong  and 
Neumann  (1982)  found  that  in  real  cases  when  the  separation  distances  are  less 
than  11  degree  of  latitude,  cyclonic  rotation  predominates.  Therefore,  the 
two  fixed  grid  points  for  the  stationary  heating  are  set  ten  degree  longitude  apart 
in  all  experiments  to  ensure  the  occurrence  of  Interaction.  After  the  dynamic 
initialization  period,  the  heating  patterns  are  allowed  to  follow  the  low 
pressure  centers.  In  Exps.  1-3  we  simulate  the  Fujiwhara  effects  in  zero 
large  scale  winds  on  a  constant-f  plane  for  strong-strong  (Exp.  1)  weak-weak 
(Exp.  2)  and  strong-weak  (Exp.  3)  storm  pairs.  Exp.  1  and  3  are  repeated  in 
Exp.  5  and  6  on  a  real  variation  of  f.  There  is  only  one  single  tropical 
cyclone  in  Exp.  4  to  help  isolate  the  effect  of  the  beta-drift.  Unlike  on  a 
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constant-f  plane  where  geophysical  orientation  is  not  meaningful,  the  inter¬ 
actions  with  real  f  for  a  weak  (west)-strong  (east)  pair  and  strong  (west)- 
weak  (east)  pair  are  quite  different,  as  we  shall  see  later,  thus  Exp.  7  is  con¬ 
ducted  to  study  the  latter  situation. 

c.  Results  on  a  Constant-f  Plane 


Fig.  7  shows  the  surface  pressure  field  at  24,  48,  72,  and  96h  for  Exp.  1. 

The  southward  displacement  of  storm  A  (west)  and  the  northward  displacement  of 
storm  B  (east)  at  24h  Indicate  that  their  interaction  has  already  caused  the  two  storms 
to  begin  to  rotate  cydonically  in  spite  of  the  stationary  heating.  The  merging 
of  the  two  storms  progresses  with  the  merging  of  the  outer  isobars  as  observed 
(Fig.  1).  By  96h,  only  the  996  mb  isobars  show  two  separate  low  pressure  centers. 

The  pressure  at  the  center  point  of  the  model  decreases  by  10  mb  while  the  approaching 
of  the  two  initlaj.  low  pressure  centers  between  24-96h  can  only  account  for  a 
pressure  decrease  of  2  mb.  This  Indicates  that  the  mutual  rotation  and  merging 
involve  dynamics  more  complicated  than  merely  advective  processes. 

Exp .  1-3  are  integrated  with  a  constant  f .  Therefore  the  orientation  has 
little  meaning  and  the  results  are  independent  of  the  absolute  initial  positions 
of  the  storms.  Fig.  8  shows  the  trajectories  of  the  storm  centers  in  Exp.  1, 
in  which  two  strong  model  tropical  cyclones  are  of  the  same  strength.  The 
trajectories  show  that  the  two  storms  rotate  about  each  other  in  a  cyclonic 
fashion  before  the  coalescence  at  102h.  The  two  trajectories  are  symmetric 
about  the  center  of  mass,  which  coincides  with  the  center  of  the  model  domain. 
Superimposed  on  the  symmetric  rotation  is  a  convergence  of  the  two  tropical 
cyclones.  The  distance  between  the  two  storms  decreases  from  ^1024  tan  at 
24h  to  ^612  tan  at  96h.  The  symmetry  remains  until  102h  when  the  two  heating 


patterns  overlap  and  one  single  large  area  of  low  pressure  is  formed. 

Exp.  2  is  identical  to  Exp.  1  except  that  the  heating  rate  is  reduced 
by  one  half.  The  cyclonic  trajectories  (Fig.  9)  are  still  remarkably  symmetric 
about  the  center  of  mass.  Because  of  the  weaker  heating,  two  identif icable 
centers  still  exist  at  120h  when  they  are  only  MOO  1cm  apart.  We  note  again 
that  at  small  separation  distances  feedbacks  between  cumulus  convection  and 
the  strom  pair's  interaction  have  been  masked  by  the  prescribed  heating  in 
our  model. 

The  speeds  at  which  the  two  tropical  cyclones  in  Exp.  1  and  2  rotate 
around  and  approach  each  other  are  shown  in  Fig. 10.  The  tangential  velocity 
of  the  cyclonic  rotation  in  Exp.  1  increases  from  Mm  s~*  at  24-36h  to  well 
over  6  m  s-^  after  72h  as  the  separation  between  the  two  tropical  cyclones 
becomes  small.  The  rotation  speeds  in  Exp.  2  are  about  1  m  s~^  slower  than 
those  in  Exp.  1.  However,  the  rate  of  convergence  seems  quite  independent  of 
the  combined  strength  as  indicated  by  the  radial  velocities  in  Fig.  10.  The 
faster  rotations  between  stronger  pairs  are  evident  observationally  at  small 
separation  distances  (Dong  and  Neumann,  1982).  At  larger  separation  distances, 
this  relationship  is  not  clear  because  the  observational  data  contains  environ¬ 
mental  Influences. 

Fig.  11  shows  the  trajectories  of  the  two  storm  center  in  Exp.  3,  in 
which  the  maximum  heating  rate  is  storm  A  is  only  half  that  of  storm  B.  We 
see  that  the  two  storms  still  rotate  about  each  other  cyclonically  and  that  they 
still  move  toward  each  other.  However,  the  trajectories  are  asymmetric  and 
the  weaker  storm  A  moves  much  faster  than  the  stronger  storm  B  in  a  way  similar 
to  that  of  a  binary  celestial  system  in  which  the  two  bodies  have  different  masses 
The  "mass''  of  a  vortex  is  perhaps  best  expressed  as  the  product  of  its  mean 
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angular  velocity  and  the  square  of  an  effective  radius  u)R  .  We  will  discuss 
this  further  later.  This  type  of  interaction  between  storms  of  different  in¬ 
tensities  has  been  observed  between  Typhoons  Flossie  and  Grace  in  1950  (Liu 
and  Wang,  1966).  Instead  of  being  stationary,  the  center  of  rotation  moved 
within  a  small  area  defined  by  the  lines  connecting  the  two  storm  centers  at 
different  times  as  in  Fig.  10. 

Results  in  Exps.  1-3  demonstrate  that  on  an  f-plane  with  no  large  scale 
wind,  interactions  between  two  tropical  cyclones  cause  the  two  storms  to  ro¬ 
tate  cydonically,  to  attract  each  other  and  to  coalesce  eventually. 

In  order  to  examine  the  momentum  fields  associated  with  the  interaction 
we  transform  the  model  500  mb  wind  fields  in  Exp.  1  onto  a  polar  grid  with 
respect  to  the  center  of  the  model  domain.  We  now  define  the  azimuthal  mean 
velocity  as 


v  )  d9 
z 
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where  v  ,  v  ,  v  are  radial,  tangential  and  vertical  velocities  on  the  polar 
r  o  z 

grid,  and  6  is  the  azimuthal  angle.  Fig.  12  shows  the  mean  vertical  (upper), 
tangential  (middle),  and  radial  (lower)  velocities  for  Exp.  1  at  24,  48,  72, 
and  96h.  It  is  Interesting  that  the  mean  momentum  fields  relative  to  the  center 
of  domain  shown  here  are  similar  to  those  in  weak  but  intensifying  tropical 
disturbances  (Hawkins  and  Rubsam,  1968).  For  example,  at  24h  there  is  a  max¬ 
imum  mean  tangential  velocity  of  ^  m  s'*  at  r  ^600  km  and  a  minimum  of 
'v-2  m  s"1  at  r  ^400  km,  reflecting  the  cyclonic  wind  fields  about  the  two 
storm  centers.  The  maximum  tangential  velocity  gradually  increases  to  M6  n  s'1 


and  moves  toward  the  center  to  a  radius  ^350  km  at  96h.  The  maximum  inflow  also 


develops  from  2  n  s'1  at  r  *^600  km  at  24h  to'U  m  s“l  at  r  ^350  km  at  96h. 
The  evolution  of  the  mean  vertical  velocity  includes  the  increasing  magni¬ 
tude  and  contracting  radii  of  maximum  upward  motion.  The  development  of  the 
mean  circulation  is  accompanied  with  the  pressure  decrease  of  10  mb  at  the 
domain  midpoint  from  24-96h,  as  discussed  before. 

The  development  of  the  azimuthal  mean  circulation  can  also  be  illus¬ 
trated  by  comparing  the  kinetic  energy  (KE)  of  the  mean  velocity  (KEM)  and 
the  KE  of  the  eddy  velocity  (KEE),  where 
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As  shown  by  Fig.  13,  the  KEE,  which  can  mostly  be  attributed  to  the  circu¬ 
lations  around  the  two  centers,  reaches  a  quasi-steady  state  after  36h.  Mean- 
while  the  KEM,  representing  the  strength  of  the  mean  circulation  as  depicted 
in  Fig.  10  around  the  center  of  rotation,  steadily  increases  until  the  coales 
cence  of  the  two  tropical  cyclones.  The  ratio  KEE/KErt  decreases  from  V5  at 
24h  to  less  than  1  at  96h. 


These  analyses  suggest  that  a  mean  circulation  relative  to  the  center  of 
rotation  develops  due  to  the  interaction  of  two  tropical  cyclones.  This  mean 


circulation  includes  tangential,  radial  and  vertical  components  resembling 
those  associated  with  tropical  cyclones.  It  is  therefore  not  surprising  that 
the  trajectories  of  the  two  interacting  storms  are  similar  to  the  trajectories 
in  the  hurricane  boundary  layer  (e.g.  Anthes  1982,  Fig.  4.6).  Compared  with 
the  nondivergent ,  barotropic  experiments  in  Section  2,  it  seems  that  the  diabatic 
heating  in  two  storms  plays  a  cruical  role  for  the  merging  of  the  two  storms. 

In  additional  numerical  experiments,  the  surface  friction  was  suppressed 
to  test  the  frictional  effects  in  the  interaction.  Results  from  these  experiments 
were  nearly  identical  to  those  presented  for  Exps.  1-3.  We  also  halved  the  coef¬ 
ficients  for  the  internal  ^dissipation.  For  the  same  'given  heating  rates,  the 
interactions  are  nearly  the  same  except  for  faster  rotation  rates,  because  the 
model  cyclones  were  stronger  with  less  internal  friction.  Therefore,  neither  sur¬ 
face  nor  Internal  friction  seem  to  be  criticAl  processes  in  the  Interactions. 

d.  Relation  of  Mutual  Rotation  Rate  to  Bulk  Parameters  of  the  System 

The  results  presented  so  far  indicate  generally  that  the  rate  of  the 
mutual  cyclonic  rotation  depends  cm  the  strength  of  the  binary  system  and  the 
separation  distance  of  the  two  interacting  storms.  Perhaps  by  relating  some 
extent  parameters  through  laws  governing  solid  body  rotation,  a  simple  description 
of  the  numerical  results  is  attainable. 

We  now  consider  the  rotation  of  the  binary  storms  is  similar  to  that  of 
a  dumbbell.  The  equation  of  motion  for  rotation  states  that  the  torque  T  acting 
on  the  binary  system  is  equal  to  the  product  of  the  rotational  inertia  of  the 
system  I  and  the  angular  acceleration  <L  with  respect  to  the  axis  of  rotation,  i.e, 

T  -  I/o  (12) 


If  we  let  R  be  an  effective  radius  of  the  storm  and  H  be  the  scale  height,  the 


mass  of  one  storm  can  be  approximated  by  Pqt:R2h.  Because  the  radius  of 
the  mutual  rotation  is  about  L/2,  the  rotational  inertia 


1  -  PQ  (K“  +  Rg)HL2  (13) 

where  Pq  is  a  reference  density  and  subscript  A  and  B  are  pertinent  for  storms 
A  and  B,  respectively. 

The  torque  is  equal  to  the  cross  product  of  a  force  and  radius  of  rotation. 

The  force  involved  in  the  interaction  can  be  approximated  by  the  advection,  then 

it  can  be  scaled  by  pQ  Y2/l,  where  2  is  the  velocity  of  the  mean  circulation 

defined  by  (8).  Because  the  mean  circulation  depends  on  the  combined 

strength  of  the  two  interacting  cyclones,  therefore  V  «  v  +  v  ,  where  v  and  v 

—  A  'B  A  B 

are  the  mean  wind  speed  within  the  effective  radius  R^  and  R^,  respectively. 
Thus,  the  torque  is  proportional  to 

T  « LH  p  (R2  +  (v  +  v  )2  L"1  (14) 

0  A  B  A  B 

We  note  that  the  torque  contains  the  dimension  of  the  kinetic  energy  of  the 
two  storms,  which  is  ultimately  related  to  the  applied  heating  Q  in  our  model. 
Substituting  (13)  and  (14)  into  (12)  and  dropping  the  over-bars,  we  get 

(vA  +  v^)2  *  (i>L2  ^  fuL2  (15) 

In  above,  f~*  is  selected  as  the  time  scale,  so  that  u  ^  ui/T  ^  fw. 

For  our  purposes  of  examining  numerical  results  of  a  limited  domain  model 

away  from  the  equator  where  f  remain  nearly  a  constant,  the  selection  of  f“* 
as  a  time  scale  is  justifiable. 
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Relationship  (16)  states  simply  that  the 
rate  of  the  mutual  rotation  is  proportional  to  the  combined  kinetic  energy 


of  the  two  interacting  tropical  cyclone  and  inversly  proportional  to  the  square 

of  the  separation  distance.  In  addition,  the  displacement  of  one  storm  should 

be  inversely  proportional  to  its  size,  because  the  radius  of  rotation  is  in- 

versely  proportional  to  the  mass  np0R  H. 

In  applying  (16)  to  the  numerical  results,  the  mean  wind  speed  within  the 

radius  of  the  gale  force  wind  (17  m  s“*)  is  used  as  v,  and  v  .  Excluding  data 

A  B 

when  the  separation  distance  is  smaller  than  the  sum  of  two  radii  of  the  gale 

2 

force,  rotation  rates  of  binary  system  at  every  6h  were  compared  with  (v  +  v  )  / 

A  B 

2 

fL  for  Exps  1,2  and  3  (Fig.  1A).  It  is  clear  the  (16)  is  a  good  description 

2  ? 

of  the  results  Exp.  1-3.  The  rotation  rates  u>  and  quantities  (v  +  v  )  /fL 

A  b 

have  a  correlation  coefficient  of  0.81.  Therefore,  our  numerical  results  can 
to  same  extent  be  represented  by  surprisingly  simple  relationship  (16). 

It  should  be  noted,  however,  that  (16)  is  arrived  through  several  simpli¬ 
fying  assumptions.  These  include  approximating  the  mutual  rotation  of  two  vor¬ 
tex  in  the  atmosphere  by  using  solid  body  mechanics  and  excluding  the  merging  from 
consideration.  While  (16)  yields  good  correlation,  it  is  only  an  approximation 
the  rotation  component  of  the  interaction. 


«•  The  Effects  of  Variation  of  the  Coriolis  Parameter 

Exps.  4-7  were  carried  out  with  variable  Coriolis  parameter,  which  can 
produce  northwestward  drifts  of  tropical  cyclones  in  the  northern  hemisphere 
(Adem,  1956;  Anthes  and  Hoke,  1975;  Madala  and  Piacsek,  1975).  The  velocity  of 
the  drift  depends  on  the  latitude  and  the  cyclone's  circulation.  To  examine 
the  free  drift  of  a  single  tropical  cyclone  in  our  model,  we  carried  out  Exp.  4. 
As  shown  by  trajectory  C  in  Fig.  15,  the  model  tropical  cyclone  has  an  initial 
northward  movement,  but  changes  toward  the  northwest  after  36h,  similar  to  the 
results  of  Anthes  and  Hoke  (1975).  The  0-72h  mean  drift  velocity  is  1.18  m  s-1 


toward  the  west  and  1.37  m  s~  toward  the  north.  The  center  at  72h  is  ^6 
to  the  west  and  ^'6. 5°  to  the  north  of  the  initial  position. 

The  latitudinal  variation  of  the  Coriolis  parameter  has  a  pronounced 
effect  on  the  trajectories  of  the  two  interacting  tropical  cyclones.  The  tra¬ 
jectories  of  the  two  tropical  cyclones  with  equal  strength  in  Exp.  5  are 
shown  in  Fig.  15.  The  two  storms  merge  much  faster  than  Exp.  1  due  to  the 
faster  northwest  drift  of  the  storm  located  to  the  south.  At  87h  only  one 
large  low  pressure  center  is  identifiable,  while  in  Exp.  1,  two  low  pressure 
centers  still  existed  at  96h  (Fig.  7).  Instead  of  rotating  around  the  + 
point  in  Fig.  8  as  in  Exp.  1,  storm  A  moves  toward  the  southeast  then  quickly 
turns  toward  the  northeast,  while  storm  B  rapidly  moves  northwestward  and  ro¬ 
tates  cyclonically  with  respect  to  storm  A.  The  two  storms  eventually  merge 
into  one  at  87h,  with  storm  B  having  traveled  a  much  larger  distance  from  its 
initial  position  than  storm  A.  The  relative  trajectories  of  A  and  B  with  res¬ 
pect  to  trajectory  C  (Exp.  4)  are  computed.  The  resultant  relative  trajectories 
(not  shown)  are  nearly  the  same  as  those  in  Exp.  1  (Fig.  3),  indicating  that 
the  trajectories  A  and  B  in  Fig.  15  is  nearly  a  linear  combination  of  the  trajec¬ 
tories  in  Fig.  3  and  the  beta  drift. 

Exp.  6  is  to  be  compared  with  Exp.  3,  where  storm  A  is  weaker  than  storm  B. 
The  trajectories  of  the  storm  centers  in  Exp.  6  (Fig,  16)  again  appear  very 
different  from  those  in  Exp.  3.  The  stronger  storm  B  shows  more  noticeable 
northwest  drift  than  in  Exp.  3.  The  weaker  storm  A  rotates  cyclonically  toward 
the  southeast  at  a  much  reduced  rate  and  with  a  smaller  radius,  apprently  due 
to  the  counteracting  beta  drift. 

Most  interesting  is  Exp.  7,  in  which  storm  A  is  stronger  than  storm  B. 

From  24  to  72h,  the  weaker  storm  B  moves  cyclonically  relative  to  storm  A. 


In  the  meantime  storm  A  moves  slowly  toward  the  southwest  nearly  perpendicular 
to  and  away  from  storm  B.  The  trajectories  take  a  strange  turn  after  72h  be¬ 
cause  the  two  storms  now  are  close  to  the  boundaries  and  start  to  influence  each 
other  through  the  east-west  boundaries  because  of  the  cyclic  boundary  conditions 
there. 

The  distictively  different  behavior  between  Exps .  5-7  can  be  explained  by 
examining  schematically  the  vectors  of  forces  upon  each  storm.  We  let  the 
northwest  drift  be  proportional  to  the  storm's  intensity  and  size  (Rossby ,  1948; 
Adem,  1956)  and  the  force  due  to  the  interaction  be  proportional  to  the  combined 
strength  of  the  binary  system  but  inversely  proportional  to  the  strength  of  indi¬ 
vidual  storm  as  discussed  in  Section  3c.  Figure  17  shows  the  vectors  and  the 
resultant  directions  of  movements  for  storm  A  at  24h  of  Exps.  5-7.  In  Exp.  5 
both  the  beta  drift  and  interaction  (both  the  rotation  and  convergence  are  counted 
for)  are  strong,  the  movement  of  the  storm  is  mostly  due  south  as  evident  in 
Fig.  15.  In  Exp.  6,  the  beta  drift  is  weaker  but  the  interaction  is  the  strongest 
the  movement  is  nearly  along  the  vector  of  the  interaction.  In  Exp.  7,  the  beta 
effect  is  strong  while  the  interaction  is  weak,  results  in  a  slow  movtt«ent  of 
the  strom  away  from  storm  B. 

4.  SUMMARY 

The  interactions  between  two  mesoscale  cyclonic  vortices  in  the  absence 
of  large  scale  winds  have  been  investigated  with  a  nondivergent ,  barotropic  model 
and  a  three-dimensional  model.  Model  results  indicate  that  the  interactions 
between  a  nondivergent  barotropic  vortex  pair  are  very  different  from  those  ob¬ 
served  between  a  tropical  cyclone  pair,  and  that  our  three-dimensional  simula¬ 
tions  agree  with  the  observed  Fujiwhara  phenomenon. 


Two  types  of  vortex  pairs  with  various  initial  separation  distances 
have  been  tested  with  the  nondivergent,  barotropic  model.  No  mutual  attraction 
is  found  in  any  of  the  cases  tested.  The  curvature  of  the  mutually  induced  ro¬ 
tation  depends  on  the  radial  profile  of  swirl  winds  (or  vorticities)  of  each 

vortex,  the  speed  of  mutually  induced  motion  is  a  function  of  separation  -  the 

closer  the  two  vortices,  the  faster  they  move.  This  is  quite  understandable, 
because  in  such  a  nondivergent  barotropic  model  the  two  vortices  can  only  inter¬ 
act  by  advection.  These  numerical  experiments  suggest  that  the  observed  Fujiwhara 
phenomenon  is  caused  by  a  more  complex  mechanism  than  just  vorticity  advection. 

Our  simulations  with  a  three-dimensional  model  reproduce  observed  Fujiwhara 
effects.  The  trajectories  of  simulated  strong-strong,  weak-weak,  and  weak-strong 
tropical  cyclone  -pairs  on  a  constant-f  plane  all  consist  of  cyclonic  rotations 
and  mutual  attractions.  The  rotation  rate  between  two  strong  tropical  cyclones 
is  generally  faster  than  that  between  a  weak  pair.  The  rate  of  convergence  of 
a  weak  pair  is  not  slower  than  that  between  a  strong  pair. 

Additional  analyses  show  that  as  the  tropical  cyclone  pair  start  to  interact, 
there  forms  a  mean  circulation  about  the  center  of  mass  of  the  two  storms  as 
the  pressure  there  decreases  more  than  can  be  expected  by  simple  advective  merging. 
The  development  of  the  mean  circulation,  consisting  of  a  cyclonic  tangential  flow 
and  a  inward  radial  flow,  resembles  the  circulation  in  weak  but  intensifying  trop¬ 
ical  disturbances.  The  kinetic  energy  of  this  mean  circulation  grows  by  a  factor  of 
four  in  72h  in  one  experiment,  while  the  kinetic  energy  of  the  circulations  asso¬ 
ciated  with  individual  tropical  cyclone  remains  relatively  unchanged.  It  suggests 
that  the  development  of  a  mean  circulation  on  the  vertical-radial  plane  relative 
to  the  center  of  mass  of  the  interacting  storm  pairs  is  crucial  in  generating 
the  cyclonic  mutual  rotation  and  merging. 

A  simple  analysis  points  out  that  the  displacement  of  one  tropical 

cyclone  interacting  with  another  is  proportional  to  the  combined  strength  of  the 


vortex  pair  and  inversely  proportional  to  its  own  size  and  to  the  square 
of  the  separation  distance.  Our  model  results  fit  this  description  well  ex¬ 
cept  for  cases  when  interacting  storms  become  highly  asymmetric  about  their 
own  centers. 

The  latitudinal  variation  of  the  Coriolis  parameter  (beta  effect)  has 
a  large  influence  on  the  trajectories  of  the  interacting  storm  pairs.  The 
beta  effect  causes  a  northwest  shift  and  a  faster  merging  of  the  two  tropical 
cyclones  of  equal  strength.  The  trajectories  of  two  interacting  tropical  cyclones 
of  equal  strength  have  a  northwest  drift  superposed  on  the  symmetrical  tra¬ 
jectories  found  on  the  constant-f  plane.  Observation  studies  showed  that 
typhoon  pairs  sometimes  drifted  away  from  each  other  if  there  were  strong  shears 
in  large  scale  flow  (Liu  and  Wang,  1966;  Dong  and  Neumann,  1982).  This  study 
Indicated  that  differential  beta  drifts  can  also  cause  the  two  interacting  tropical 
cyclones  of  different  strength  to  diverge  when  the  one  initially  located  to  the 
west  is  stronger. 

These  findings  should  not  be  accepted  without  caution  because  of  several 
limitations  of  the  numerical  model.  The  model  domain  is  perhaps  too  small  for 
two  tropical  cyclones.  In  addition,  the  horizontal  resolution  of  V3  is  only 
marginal  for  resolving  realistically  the  smaller  scale  dynamics  near  the  center. 
Being  a  uniform  grid  model,  without  decreasing  the  horizontal  resolution,  the 
model  domain  cannot  be  expanded  due  to  limited  computing  resources.  The  cyclic 
boundary  conditions  created  problems  (as  evident  in  Exp.  7)  when  two  storms  may 
have  interacted  with  each  other  through  the  east-west  boundaries.  Perhaps  the 
most  serious  limitation  of  our  simulation  is  the  heating  prescribed  ji  prior  in 
the  three-dimensional  simulations,  which  may  have  masked  the  interactions  between 
two  adjacent  tropical  cyclones  on  the  scale  of  cumulus  convections.  However,  the 
development  of  the  mean  circulation  about  the  center  of  mass  of  the  two  tropical 


cyclones  occurs  at  very  early  stage  of  the  interaction  when  the  separation  is 
still  large.  This  suggests  that  the  detailed  characteristics  of  cumulus  con¬ 
vection  in  individual  storm  may  not  be  important  in  setting  up  the  cyclic  ro¬ 
tation  and  mutual  attraction.  The  use  of  the  prescribed  heating  was  justifiable 
except  at  small  separations  where  the  divergent-convergent  pattern  in  each  storm 
s>ay  be  modified  due  to  the  proximity  of  another  one. 

In  future  research,  a  parameterized  convective  heating  should  be  utilized 
to  investigate  the  abovementioned  secondary  effect  of  the  cumulus  convection. 

In  addition,  the  parameterized  heating  may  react  to  large  scale  winds  in  a  non¬ 
linear  fashion.  Therefore,  the  nonlinear  effects  of  the  large  scale  winds  on  the 
Interactions  of  two  tropical  cyclones  also  ought  to  be  studied.  The  question  of 
what  is  the  maximum  separation  distance  for  storm  pair  to  interact  is  also  left 
for  future  studies  when  numerical  models  of  tropical  cyclone  cover  a  larger 
domain  are  constructed. 
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FIGURES 


Figure  1:  Surface  isobaric  analyses  at  ooooZ  on  Sept.  15-20,  1964  showing  the 

rotation  and  merging  of  typhoons  Kathy  (K)  and  Marie  (M).  Isobars  are  plot¬ 
ted  every  3  mb.  (From  Liu  and  Wang,  1966). 

Figure  2:  The  radial  distribution  of  relative  vorticity  for  type  A  (solid)  and 
type  B  vortex  (dashed) . 

Figure  3:  The  trajectories  of  type  A  (solid  lines)  and  type  B  (dashed  lines) 
vortices  at  separation  distances  of  300,  400,  600,  and  1000  km.  The  cross 
is  the  center  of  the  domain.  Time  interval  between  two  adjacent  dots  is  12h 
Squares  denote  the  initial  vortex  centers. 

Figure  4:  The  vertical  distribution  of  heating  used  in  the  model  plotted  on 
an  arbitrary  scale  (solid),  as  compared  with  (Tc-T)  for  mean  hurricane 
sounding  produced  by  a  one-dimensional  cloud  model  (Anthes,  1977,  Fig. 4). 

Figure  5:  The  horizontal  distribution  of  heating  used  in  the  model  plotted 

on  an  arbitrary  scale  (solid),  as  compared  with  the  observed  radial  dis¬ 
tribution  of  rainfall  rate  in  a  typhoon  (Adler  and  Rodgers,  1977,  Fig.  5  ). 

Figure  6:.  The  radial  distribution  of  the  quasi-steady  wind  speeds  in  model 
layer  six  (o  ■  0.85)  and  seven  (5  *  0.963)  generated  by  the  stationary 
prescribed  heating. 

Figure  7:  The  surface  pressure  field  at  24,  48,  72,  84,  and  96h.  Contour 

Intervals  are  4  mb,  the  outmost  closed  isobars  are  1008  mb.  Longtltudes 
are  arbitrarily  set. 

Figure  8:  The  trajectories  of  storm  centers  in  Exp.  1.  Numbers  on  the  curve 
denote  times  in  hour. 

Figure  9:  As  in  Figure  8,  except  for  Exp.  2. 

Figure  10:  The  tangential  and  radial  speeds  of  the  two  Interacting  tropical 
cyclones  relative  to  their  centers  of  mass  in  Exps.  1  and  2. 

Figure  11:  As  in  Fig.  8,  except  for  Exp.  3. 

Figure  12:  The  azimuthal  mean  radial  (vr),  tangential  (vg),  and  vertical 
velocity  (w)  of  the  wind  fields  relatives  to  the  center  of  mass  in 
Exp.  1  at  24,  48,  72  and  96h. 

Figure  13:  The  development  of  the  kinetic  energy  of  the  "mean"  flow  (solid) 
relative  to  the  center  of  mass  and  the  kinetic  energy  of  the  "eddy" 
associated  with  the  two  storm  centers  (dashed)  in  Exp.  1. 

2  2 

Figure  14:  The  rotational  rates  u  compared  with  (v^  +  v_)  /fL  . 


Figure  15:  The  trajectories  of  the  free  drifting  storm  in  Exp.  4  (Curve  C) 
and  of  the  two  interacting  storms  in  Exp.  5  (Curves  a  and  B) . 

Figure  16:  As  in  Fig.  15  except  for  Exps.  6  (solid  lines  with  dots)  and  7 
(solid  lines  with  squares). 

Figure  17:  Vectors  showing  schematically  the  force  of  interaction  (I),  Beta 
drift  (D),  and  the  resultant  movement  (M)  for  tropical  cyclones  A  at 
24h  of  Exps.  5,  6,  and  7. 
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AN  AXISYMMETRIC,  NUMERICAL  MODEL  FOR 
A  NON-HYDROSTATIC  BOUSSINESQ  OCEAN 


1.  GOVERNING  EQUATIONS 

The  governing  equations  of  the  axisymmetric,  non¬ 
hydrostatic,  Boussinesq  ocean  model  are 
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where  7“  =  — j  *  —  other  symbols  are  listed  in  Appendix  A 


Above,  the  density  anomaly  b  is  defined  according  to 


p  -  o  (=) 

b  -  - TT1 -  §  . 
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where  c  (=)  is  a  reference  density  and  is  a  function  cf 
deDth  onlv.  Brunt-Yaisalla  frequency  \'_  is  defined  as 
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The  continuity  equation  is  that  of  t'r 
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2.  THE  MODEL  GRID 

It  is  determined  that  a  fully  staggered  grid  is  most 

...  * 
expedient  tor  storage  economy  for  a  given  spatial  resolution. 

As  shown  in  Fig.  1,  the  radial  (u)  and  the  tangential  (v) 

velocities  are  defined  at  cross  points,  vertical  velocities 

(w)  are  defined  at  open  circle  points,  and  the  pressures  (PJ 

and  density  anomalies  (b)  are  defined  at  blackened  dot  point 

This  grid  system  has  the  following  advantages: 

a)  it  saves  storage  for  a  given  spatial  resolution 

b)  it  is  very  economical  in  terms  of  number  of 
computational  operations  for  the  finite 
difference  (FD)  equations  of  (1-1)  to  (1-4). 

c)  it  is  very  easy  to  specify  the  boundary 
conditions , 

d)  the  pressure  diagnostic  equation,  of  the 
elliptic  type,  can  be  reduced  to  the  standard 
form,  and 

e)  there  is  no  spatial  separation  of  solutions 
on  the  grid. 

In  order  to  consistently  index  the  grid  points,  we 
let  index  pair  (ij)  represent  the  i-th  point  in  the 
r-direction  and  i-th  point  in  the  r-direction.  In  addition, 
m  is  the  maximum  number  of  points  in  the  r- direction ,  and  n, 
the  maximum  number  of  points  in  the  c-direction.  ~he re  f ore 
there  are  m  x  (n-1)  points  for  radial  and  tangential 
velocities,  (m-1  x  r.  points  for  vertical  velocities,  an_ 
i.m-lj  x  i.n-1;  points  for  mass  distribution  <"c  and  p  . 
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Fig.  I  The  fully  staggered  grid  system  of  Lhe  ocean  model 
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4.  DERIVATION  OF  THE  DIAGNOSTIC  EQUATION  FOR  PRESSURE 

The  nonhydrostatic  pressure  at  tine  t  is  needed  t 
compute  the  pressure  gradient  forces  in  '3-2)  ana  1 3 - 5 ■ . 
To  "recover"  the  pressure  from  the  motion  fields,  we  make 
use  of  the  continuity  equation  by  differentiating  ll-". 
with  time  we  get 

—  1  +  _9_  -  n 

r  3r  r  at  32  at  “  * 

which  can  be  written  in  finite  difference  form  for  a  mass 
point  ii  as 
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Substituting  (3-2),  (5-5)  and  (4-2)  into  (4-1),  we  have 
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We  obtain  the  standard  form  of  an  elliptic  equation  in  FD  fo 

AX  i  Pi  - 1  i  *  A:jPii-i  ■"  BBij  Pii  *  txiPi*ii  *  c:jpii*i  * 

Equation  (4-5)  can  be  solved  numerically  by  the  SEVP  solver 
(Madala,  1 9 “ 8 )  ,  providing  the  boundary  conditions  are 
properly  posed. 

The  conditions  for  the  four  boundaries  are  determine 
according  to  the  following  assumptions: 

(a.  At  (r,).  =  (r,i.  *  0,  the  natural  condition  for  the 

cylindrical  coordinates  calls  for  u  =  v  =  0  =  Fu/~t  =  ?v  it, 
the  gradient  balance  requires  that  fi?  ir' r.^  =  0.  Therefor 
an  extra  column  of  ?  is  needed 
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Note  that  if  is  set  equal  to  zero,  (4-7)  describes  a 

non-divergent  and  zero-vorticity  boundary  condition  at 

r  =  (r, )  .  Once  v  .  is  determined,  a  gradient  balance  at 
v  I'm  mi  ° 

r  =  (rl^m  rec*uires 
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where  a  column  of  dummv  points  p  .  has  been  introduced  for 

mi 

computational  purposes.  The  second  part  of  the  RHS  of  ,4-S 
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At  the  bottom, 


ese 
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into  the  continuity  equation  ,'i-l  >  ,  we  get 


Following  the  same  deduction  between  (4-1)  and  (4-5.)  ,  we  get 

an  expression  similar  to  (4-5)  with  the  second  team  on  the 

LHS  and  G^.  in  the  RHS  absent.  Thus,  can  be  obtained 

by  the  same  SEVP  solver  by  setting  =  0  and  G ^  =  0. 

(d)  At  top  wwn=3/3t  *  0.  Following  the  same  line 

of  reasoning  as  in  (c) ,  we  obtain  P^n  by  solving  (4-5)  with 

C-n  pin*i  *  0  and  3  °- 

n  m+l  m 

In  summary,  the  elliptic  pressure  diagnostic  equa¬ 
tion  (4-4)  is  to  be  solved  with  the  following  boundary 
conditions 

1.)  At  r  -  0  Poj  -  Pjj  i.*..  (4-6) 

2)  At  r  *  irj.  P  .  *  P„  ,  •  *  function  (v  • )  (4-S) 


LIST  OF  SYMBOLS 


an  array  of  constants,  varying  only  in  redirection, 
defined  by  (4-4),  used  in  (4-5) 

an  array  of  constants,  varying  only  in  e-direction 
defined  by  (4-4),  used  in  (4-5) 

an  array  of  constants  related  to  r^  and  Ar^  used  in 
(4-2) 

an  array  of  constants,  used  in  (4-5) 
density  anomalies,  defined  in  (1-5),  cm  s  ^ 

’  an  array  of  constants,  varying  only  in  r-direction, 
defined  by  (4-4),  used  in  (4-5) 

an  array  of  constants,  varying  only  in  z-direction, 
defined  by  (4-4),  used  in  (4-5) 

an  array  of  constants,  related  to  r^  and  Ar^  ,  used 
in  (4-5) 

Coriolis  parameter,  s  * 

_  ~> 

gravitational  acceleration,  cm  s 

an  index,  denoting  i-th  point  in  r-direction 

an  index,  denoting  i-th  point  in  z-direction 

?  - 1 

horizontal  diffusion  coefficient,  cm"  s 

•>  -  x 

vertical  ditfusion  coefiicient,  cm"  s 
left  hand  side 

the  maximum  number  of  grid  points  in  r-directicn, 
upper  bound  of  i 
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m- 1 

m  m-  2 

\'_  Brunt- Vaisalla  frequency,  s  1 

r.  the  maximum  number  cf  grid  points  in  c-direction, 

upper  bound  of  i 
n-^  n- 1 

n,  n-  2 

.  7 

p  pressure,  dyne  cm  “ 

RHS  right  hand  side 

r  radius,  cm 

radii  of  momentum  points,  cm 
T-,  radii  of  mass  points,  cm 

Ar^  distance  between  two  hori contally  adjacent  momentum 

points,  cm 

Ar?  distance  between  two  hori contally  adjacent  mass 

points,  cm 

SEVP  stabilized  error  vector  propagation 

t  time,  s 

At  time  interval,  s 


u 

radial 

velocity, 

cm  s 

V 

tangent 

ial  veloci 

ty,  cm  s 

-1 

w 

vertica 

1  velocity 

,  cm  s  * J' 

height 

from  ocean 

bottom, 

cm 

heights 

of  circle 

po ints  , 

cm 

hei ghts 

of  cross 

and  dot 

points,  cm 

16 


distance  between  two  vertically  adjacent  circle 
points,  cm 

distance  between  two  vertically  adjacent  cross  cr 
dot  points 
density,  g  cm  ^ 
a  constant  density,  1  g  cm  ^ 

a  reference  density,  varying  only  in  c-direction 
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APPENDIX  A  -  FORTRAN  CODE  FOR  THE  NON-HYDROSTATIC  MODEL 


A  listing  of  FORTRAN  code  of  the  ocean  model.  The 
major  functions  of  the  main  program  and  subroutines  are  as 
follows : 


OCEAN 

IXIT 

START 
j  PUTOUT 
|  MAP 
ADVECT 

DIFF 
I  PRESS 

'  MAT I XV 

J  BSM1 

I  BSM2 
\BSM3  • 

FRWRD 

BOUNDV 


main  program,  calls  all  subroutines,  manages 
job  flow,  controls  input/output. 

sets  up  independent  variables,  defines  con¬ 
stants 

defines  initial  conditions 

gets  various  fields  ready  for  output 

prints 

computes  all  inviscous  terms,  except  for  the 
pressure  gradient  forces 

computes  horizontal  and  vertical  diffusions 

solves  the  pressure  diagnostic  equations  and 
computes  the  pressure  gradient  forces,  appears 
only  in  the  non-hydrostatic  version 

inverts  matrices  J 

>  Used  in  SEVP  method 

matches  forward 

sets  outer  boundary  conditions  for  momentum 
checks  if  the  time  step  is  linearly  stable. 


CHECK 


CARO  Image  MtE  Er'IT8s<CI^E*'J  ••  v£RS2s>  05. 2*  C‘TE»  1 0/lt/il 


TI"E* I « 1 1 ; 1 1 2 1** 


«$rut  bce*ko,S8lrco»rr:,'t«se: 


2* 

3* 

«• 

5* 

*• 

Jt 

8* 

4* 

10* 

lit 

12* 

13* 

t«* 

13* 

I** 

IT* 

IS* 

1«* 

20* 

21* 

22* 

23* 
2«*  C 
23*  C 
2S*  C 
27* 
2«* 

2R* 

30* 

31*  C 
32*  C 
33*  C 
3«* 

33* 

3*» 

37*  C 
3**  C 
10*  c 
*0* 
at*  C 
*2*  C 
*3*  C 
a  «* 

*3*  C 
«*•  C 
*7*  C 

*0*  c 
*0* 

*o*  c 

51*  C 
32*  C 
33*  C 
!«■ 

35* 

5*« 

37.  C 
30*  C 
SO*  C 
*0* 

*1*  c 

*2*  C 
*3*  CC 
*«• 

*3* 

*0* 

*7.  C 
**•  C 

••*  c 

70* 

71*  C 
72*  C 
73*  C 
7«* 

75* 

70* 

77. 


program  ocean 
P4rAmET£w  H*21,k*;i 

parameter  mi*n.i  ,  ►2**‘«2,ni»n-i  ,N2«t.«2 
PARAMETER  Nf'»2***Al*Ml**‘*»l*Rl 
DIMENSION  DATA  I (NC),CATA2(NC),BATAltNC) 

C0MM#N/0t,£/VR|  (A,A1),VT1  (“,M),VZl(*t,iOr8l(»l,M).VR2(««,*in, 

1  VT2(m,m  ),v:2(M|fNj,s2(»j,Nl  ),V*3(«,M).*TJ(M#NJ), 

2  v23CM|,A),83C“l,Nl},K«t,Nn 

EOul VALENCE  (OATAl . Vfll )# (DATA2.VR2] . (QATA3.VR3) 

DATA  OATai/Nd*0,/,CATA2/NO*0,/,DATA3/N0*0./ 

COMMeN/TMR/RKMJ.REIPn.OPl  (“t),0R2(M),2i(N),22(Nl  >,D2l(Nl  ),0Z2(N) 

CAMM0N/7WR/RHR,fcMfR(M),8V2tN),Al,PMA,tNeA,*N08,CS»I#e,M«<M)#2*tN) 
C«MMPN/FfR/DELT,  RTJME,  I  T 1  *£ ,  I S  T EP ,  I  SmO  ,  I T  APE ,  T|  V 
C*lL  INDu«P 

190  PfflMAlClJ! 

BEAO(5,100mi»e 
RE AO (3. 100) ITER 
Rl*D(S,tOfl)I(H'T 
RE  An  (5, 1  99)  I$>*e 

ISTEPPO 

PlAO(b,109)ITAPl 
C*LL  INI T 

ip()TiMftEa.o)Gf  is  te 

CONTINUED  INTEGRATION  PRO"  A  NISTORT  Tape 
REAO(1)ITIMC.DATA1.0ATa2.P 

c?  10  2* 
ie  c *ll  start 
2o  *7l"EPI7lME*i*9fl. 

POINT  3LT  I N I T I au  FIELD* 


CALL  MUTOUT 
I*(ITtR.E0.9)STfP 
CO  00  ISTCPP1.ITER 


CP-PlTt  ALL  INVISCIB  TfR-S 


c*ll  »dvect 


call  I’tFF 


CALL  PRESS 


Compute  VISCOUS  TERMS 


ANC  ACC  P°ES3LR£  GRADIENT  FORCES  ’0  TENDENCIES 
DIAGNOSE  (RECOVER)  Tn£  PRESSURE  »IELO 


•aocmInS  In  TI«E 

FIRST  1 1 mE  STEP  IS  *  PR  *  A  pp  IF  START  IS  CAuLEC 

IF(ISTEP,E(J.l,A*C,ITI*E.EG,0)CCLT*0.5*eELT 
CAUL  *»«RD 

IF<IS»EP.EG.t .ANC,ITI*E.Ea,0)PELT*2.*0ELT 


CAUL  OPUNOV 


Of  F I nE  OmUnOaRT  VALUES  P»R  VEL8CITT 


Ct-EC*  IF  CELT  IS  STASLE 


C*LL  CmIC* 

ATI"t**TI”E*CELT 

ITi«e*«TI-e/3P0L. 


®  a  I  *  T  OUT  RESULTS  EVE»»  Put  STEPS 
P  (-OCIISTCP.IflT  1  .EG.  c  1C  ALL  PuTOijt 

-APE  -ISTORv  tape  Evert  jtare  STEPS 

I*f“flC(ISTE»,ITAPE).E:tC)MRPE(2):TlMEf0ATAl,0ATA2,P 

«e  c*ntinue 

STOP 

END 


008)000 

0002000 

00»300C 

OOOROOO 

000500C 

oouoooo 

003T00C 

oooeooe 

ooeoooc 

0010000 

oonooo 

0012000 

0013000 

001R00C 

0013000 

OOIOOOO 

0017000 

ooisooo 

OOIROOO 

0020000 

0021000 

0022000 

0023000 

002POO0 

0023000 

002*000 

0027000 

0028000 

0024000 

0030000 

0031000 

0032000 

ooisooo 

ooioooo 

ooSsooo 

ooSoooe 

0037000 

ooioooo 

ooioooo 

00*0000 

oo«togo 

00*2000 

*0*3000 

00**000 

00*30«0 

00**OgC 

00*7000 

oo*eooc 

00**000 

0040000 

OOSIOOC 

0092000 

OOSSOOO 

ocs*ooe 

0033000 

003*000 

oostooc 

003*000 

003*000 

00*0000 

00*1000 

00*2000 

00*3000 

O0«*00C 

00*3000 

00**000 

OOeTPOC 

oceeocc 

00*ROOC 
OOTOCCC 
CC'ICSO 
00  ?2CC0 
0073000 
0C78C0C 
007500C 
007*000 

ooTTooe 


C*«0  1**61  PJUE  £0IT8«fCl»£<*)  —  ve«SlC»  03,2*  o» re*  i  u/2* /ez 


121121** 


•  ••  •>£»»£«  IK  I  f 

l* 
i  t 

3* 


ft 

r 

8 

4 


13 

1« 

IS 

1ft 

I  7 
ia 
i' 

«o 

21 

22 

23 

2* 

2S 

2ft 

V 

a 

2* 

JO 

II 
32 
13 
3« 
JS 
3* 

37 

38 
38 
«0 
al 
«2 
a] 
ft» 

«5 

•  ft 
«7 
at 
a* 
80 

51 

52 

53 
5* 
35 
3* 
Jl 
58 
38 

to 

tl 

*2 

*3 

** 

*5 

ft* 

*’ 

*8 

ft* 

re 

71 

72 

73 

7* 

73 

7* 

77 

78 
7ft 
(0 


S'-8B8*'  T  I  k£  X**I7 

PftBftwtTE#  »«21,k*2I 

P*ft*"t7£h  “!•“•)  k i i 

C»«*-9*»/eN£/Vftt(kfKiJ,wtHK#Xt),'<2l  (*l,*},81(*l,ftt),V»2(»,»  l), 

1  v72{*‘,Ml,v2i<«M,KJ,a2C'M,Kll,V»!(*«,*l),vTJ(ft,K!), 

2  vZS(“i,K),8J(>*l,*l),B(«i,ftt) 

C'*<“8*»/T.ft/B|  («),«2(»l  l.osi  t-l  }.OR2(»),21(N),22(Kl  ),CZl  (Ki  ),0Z2C*) 
C8K»<8N/TKft/8H8,k>.f)j(M  >»8»'2(K;,fti.»Ki,eKOft,8,-o*(eefti,s,Ki(  (“ifZuc-i 
Cft«.«8K/7  8#/C£LT,*1  t**£,I  TI<.e,I9T£»,I3K0.IT*8E,T*V 
P  *#ft**t  TE#  *#•“♦!, KP*K»I 
Pft#**ETEP  NPu1'*? I  »*‘8l«»l 

P*Bft»ETEfl  , »P2ft«P.j,HPi»HB.i  ,KP2ftKP.; 

ft£it»ft  ftcstfOlKV/SisvijPTitn* 

C»»»8K/t»P/»IMVt«P2,K*J.NeL«)fftlKW| (KP2,KP2,KeLP|),»Cfftf"*/3)» 

1  «TlLf  ft(P82l,P(“P*KP),.M8$IZ2(ll<Bt«),IStKei.«)»»tP»(*'eU*>, 

2  IE  (**(.•  )«P  It  <**8)  »P  IK  C* )  #P2l(*i8).P2«  (ft#  },**{»#),  ft  »  (ft#) » 

3  8e<**,KB).C*t*«8)«Cr(*0) 

C 

C  riJTjftLlIC  ftUL  OEPCKI5EKT  VftOTftBcCS  »NC  CftKSTftftTS 

e 

c 

C  *l*K*  13  The  *e*nt#EftSIl!N*U  S^SCr-IKC  C SEP , 

C  BPS  TIP£  SPfte’HINS  IK  SIJ8»8UT1k£  PB»#0 

c 

rElTftftOO. 

ftL*»ft*0,IO 

6*880, 

L»T»IO 

C9»Ia2,ft7.2722£.s*slN(tftTft].ia!5ft/|8o. ) 

c 

C  OfPIftE  ftftOt!  ft 7  6»I0  P8IKTJ  ftKQ  «LL  6*ID  JKTE8»»LS 

C 

fP  10  I«l,Mt 
10  r«l  (I  lft2ii.E5 

tit  20  1*2, N 

20  #1(1)081  (t»l)»OM(I*t) 

0*  30  1*1, **l 

30  82(»>»o.3P(*m)«*lU«t)> 
a#au)o2.ft(#2(i)-Ki(tn 
3#2*,,)o2.ft(#l  (* )•# | (* 1 ) ) 

0*  ««  I*2<M1 
»0  u#2U)o»2(D«»2(I*l  J 
C  ***■•*•  i"  *0  ("#  1  j 

C  5*«ftioo*t  (•*»* ) 

e  oi p i#£  *ll  pz’s 

r,t  1«1  J«i 

mo  2  <  (  J  )•  (  J»1  1  *200,1  a 
•l’  t|0  Jftt,".  ( 

110  5Zt (J)a2l  (J*| 3-2  7 (J) 

CZ2< I  1*2, • (0,3*1  Zl (1 l*Z'  (l j) 

22(1)00. 5*022(1  j 
00  120  J*2,*l 

0Z2U)a0,3*(0Z!  (J)*CZl  (J.l)) 

120  Z2(J)»Z2(J-1)*022(JJ 

nZ2(Mo2,*(Zi  (K)»Zl  (K-l ). 0,3*02'  f •* |  )) 

C  «**o»i**p*6(02i ) 

C  3Z"ft»oOZl (Ofti ) 

C 

C  if  P 1 KE  C^llftftTJ  PPB  S£V#  S?LV£# 

C 

**(1)00. 

CM-PloO, 

*»(!)• 0. 

C',f<*l»0, 
s*  so  i»i,*t 

»«(!•!  >•#)  (I)/(o,ft*(#I  (I)*81(I*U)*C8I  (I)*(5#2U)} 

30  C> (1*1 )•#! (1*1 )/(C,3*(#t (t)»#l (1*1)  )*C8  1  (I J*0#2(I*I!  5 
C"  *0  J*l,->| 

»T C j*l )8T./(02l(J)ftC22(J) ) 
to  c’(J*i)*i«/(ozi(j)«rz2fJ*')) 

P*  7.1  J«1,N# 

58  70  I»i,m» 

70  M(I,J)».C*  (I  )*ftl(  I  >*C7  (.*)•*»(,’) 
sSiZEoki /K«t*  *«pc (n l ,2  j 

e*  bo  HP® i , Kftco i 

#0  K#»:Z2(S8)oKSt2E 

K8|lZ2(K9lK)o’n*(*Bl**1  J-KJIZSftKOL*  1 


o  o  c  i  c  o  c 
o  0  o  2  c  u  c 

00U3C00 

OOUftOoO 

00U5C0O 

ooo*ooc 

0007000 

0008000 

0008000 

0010000 

0011000 

0012000 

oouooo 

OOlftOOO 
0015000 
001*000 
00  17000 
00  I  8000 
0018000 
0000000 
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0002000 
0003000 
OOOftPOO 
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0007000 
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OOOftOOO 
0030000 
0031000 
00)2000 
00)3«00 
0034000 
00)3000 
OOlftOOO 
00)7000 
00)8000 
003*000 
00*00g0 
00*1000 
00*2000 
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00*ftOoO 
<10*3000 
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00*8000 

ooSocoo 

0  03 1 0  00 
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0033000 
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OOftlOoC 
OOftftOoC 
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OOftftOoC 
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00*8000 
OOftftOOO 
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OOMOoC 
00  72000 
OO73CC0 
007*000 
0073000 
007*000 
0077900 
0071000 
OO'ftCOC 
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*  ISO  fe  i»E  C"nST*nt  LSft  !*■  Subroutine  e*L*S* 
FfS  COSJTlsT  PJV.  *se  y'RT,  CONDITIONS 

?nc**n  t  €•* t ) /s  1 C*" ) 

osc*>»mt  (**i  )*®i  (»))«o*t  («t  )/((Rt  ft  )*"t  (*<2)  )*«t  (*<)*0Bi(t>2n 
fS08»0, 


CtrtSE  CENSITT  B£LiT£C  CONSENTS 


«*-o*i. 

0*  t)0  J*1,N 

11*  OVjtj j«|  ,£•* 

TBv»0. 

0*  1 33  J»1,n 

135  T8v«***XI  1T8V,8*2(J)  1 
T3v»t./S8BT(T5*J 

Cl  FINE  M»IZONT*t  *ND  VERTICIL  CIFrujIOS  J 6£FF JEN T 3 

C"EF"»0,002«DB| ( 1 J**2/0tlT 
C0EFZ«0.0OI*02|(l )**2/CElT 
0"  1«0  X«t,**| 

ISC  h« (I )sC0EFh*{1 ,»5a*tiR {•FLfllTt ««•!)/?.) ) 

DO  ISO  Jal,N| 

ISO  ZF(J)«C"EFZ*(t.«3.*(E«*<»FL9*TlJ-l)/S.)«E*F(-FL6»T(Sl»J)/5.))) 

RETURN 

ENC 
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2*«  p<*  2o  1*2, mi 

27*  2o  vT8d,J>*«0.25*t {VR tI,J)AWR(I-t»J))«(YT(I,J7.VTtlAt .J)),CPl(I»l ) 

26*  1  ♦(VRClAttJ)*YR(I,J})A(VTtI*t,J)*YT(I,J))/p*ttI)} 

2*»  2  AVTJCI.J) 

3«*  C 

3t*  C  m9HZ6NT*l  A0VECTI8N  F9R  VERTICAL  VELBCITV 

32*  C 

33*  P*  30  J*2,nt 

36*  08  30  1*2, M2 

3S*  30  VZ3U,  J)»-o.25a(  (VR(I,J)*VR(I,J*I  ))*  ( VZ C 2 , J1 -VZ Cl  *  1 «  J)  )/PP2(I) 

36*  1  *(VR(i*i,j)*vR(i*t,j.i>)<(VZti*l,J)>vza*2))/0R2M*n) 

37*  2  *vZ3(i,j) 

38*  00  «o  J«2,ni 

JR*  as  *Z3d»J>*-0.25*<vt.<2#J)*vR(2,J-1 » >•  (VZ  (2,  J)-vZ  ( 1 , J) >/DR|  (2) 

«0*  t  AVZJtl.J) 

«1*  DO  So  J«2,nt 

«2*  So  vZ3tMt,J)*-o.2S*(vR(Mi ,J)*v*(«l,J.t ))AfvZt"t,J)-vZ(“2,J))/CR2(»t) 

03*  1  «VZ3(»1,JJ 

«*•  C 

*5*  C  H9N I28MT *L  APvECTJBm  FOR  BUOYANCY 

lit*  C 

*7*  DO  NO  J*1,Nt 

«8*  CC  NO  1*2 ,m2 

«R*  80  e3(I,J)«e3(I,J)-fl,«*tVR(l,j).(9(I,J).i»(I.|,J))/OR2(IJ 

58*  1  ♦yR(I«l  ,J  J*<9(t*1 ,  J)-9(I,«’J  )/0"2(I*t  ) ) 

St*  C8  7 o  J* | ,Nt 

52*  7o  B3(t, 21*93(1, J)*8,S*vR(2,j)*(B(2,J)*B(l,J)>/DR2(2) 

53*  CO  *0  J* 1 , N t 

5«*  *0  6S("t,JJ*<*SCMl,J>*0,5*v*(>‘i,J).(B(Mt,Jj»B(M2,J))/CB2(Ml) 

55.  C 

58*  C  ,  VERTICAL  ADVECTION  FOR  RADIAL  VELBCITY 

57*  C 

58*  pe  80  J*2,N2 

58*  po  *u  1*2, Ml 

*8*  8o  vRS<lAJJ»VR3(I,J>.0.23*((VZU-l*J)*VZ(I.JI)*{V«U,J)-VR(I,J.n) 

8 1  •  i  /oz2cj>*(vzii,j*i)*vz(i»i,j*nt*(vR(i,j*i)«vR(i,jn 

•  2*  2  /nZ2( J*  1 ) ) 

*3*  C"  *5  1*2, **t 

8  4*  8*  vP3d»  1  )*VR3(I,  l  )*C.25*  (VZ  (I,2)*VZC»1 ,2>  )*(VR(I,2)«V*(I,  t  n 

85*  1  /DZ2  (2  > 

68*  o*  *e  2*2, •>! 

07*  8 1  vP3d*M)*vR3d.‘>1)*C.25*(vZ(:-i.Nt)*vZdtM))«(y8(I,sti*vRd,N2j} 

*8*  -  I  /0Z2  (M  ) 

6  8  •  C 

70*  C  VERTICAL  ACVEC7I0S  F"R  TANGENTIAL,  VEL5CITV 

7  1  •  C 
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0050000 
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??• 

73* 

7  «* 

75« 

7** 

77  • 

7I« 

79* 

*9* 

Cl* 

C2* 

8  J*  c 
««•  c 

C5*  C 
C6* 

•  7. 

ee* 

C9* 

90* 

91*  C 
92*  C 
91*  C 
9«* 

95* 

99* 

97* 

98* 

99. 
190* 
101* 
m2*  c 

103*  C 
10«*  C 
1 03* 
149* 
107* 
1«8* 

1 09*  C 

no*  c 
111*  c 
112* 

111* 

U** 

115*  C 
119*  c 
t|7*  C 
1 1 •* 
in* 

129* 

121* 

122* 


loo 

■  ■  i  uo  i •  a .  “ i 

loc  .7j[i,j)*vT5(t,j).o.25*((vz(i.i<j).v2(j,j})*(vT(i#j).WT,jj.ljj 

1  /nl2(JJ*(V2C!,J»l)*V2CI»!|J*l))*(VTn#J*l).yT(*,jn 

2  /P2«(J«m  * 

I*  107  !«2,»| 

io«  » t  j{ i, t  )»v71(l , » ).c,25* tvl(!,2i«vr (1*1,2) ) • (VT (1 ,2 )-vT ( i ,  i )) 

'  /P  Z2 ( 2 ) 

C*  10*  1*2, “1 

10*  otJtl.M  )*VTJ(I,I.|  J«.>,25*c*2(ift.i)*v2n.i#si))*cv7tl,*.|  j.vT  1 1  .*,» }) 

1  /(122  f*>  I  1 

vf»tJC*L  ooveCtija  rg*  »E*Tlc*t  VELOCITY 

C*  1*0  J*2,M 

i«o  t*i,»i 

l««  02ia#J)*v2J{t,j).n.25*ttvZU,j.i)«vzti,jj)*(v2t:,j).wz(j,j.n  ) 

'  /"2l (J-i)«<*2(I, J*l)*YZ(I,J)).CVZ(I,j*i).yZ(j,  . j } 

2  /OZHJj) 

Yl  *7  IC*L  *CvEC7!9N  *8*  B 

.7"  is*)  J»2,'<2 

£9  150  1*1, “1 

150  a!fI»J)*8Jl!,J)»o,5*£v2fI,J)*{8Cl,jJ«B(I,j«t))/(>22{j) 

1  ♦yZU,J*1)»<BCI,J»i)«B(I,J))/cz2(J»ij) 

P*  1*0  t*l,"l 

i*o  9JU*  l)*9  JU,tl.u.5*Y2(l,2)*(«2fl,2)*82(l,lJ)/0  22(2) 

I7« 

17c  *3<  1  ,'M  )«B  J(I,M  )«*,5*v2  (1,»M  )*(4  (I  ,M  )•«<  t  ,NJJ  J/CZ2  ( K 1  ) 

INC 97  Z*  Tt»N$  889  H09IZ8NTAL  n8hEn tu* 

C."  110  J9  1  ,  N  t 
C«  110  I«a,-l 

V9J{I#J)«V81(I,j)*vT(I,J)«(V7(I,J)/*1  (I)*C8»n 
tic  VT3(I(J)4VTl(t,J>.v4(I,J)«(VT(I,J)/4}(J)*C94t) 

BcCy »nCy  tt»«  Mil  yErt ICac  *CCEL£**TJ5n 

C"  120  J«2,*1 
07  120  1*1, -I 

120  Y2jti,jj*vz3(i,j).o.5*(»a,j)*ari,j.in 

STAAI 19 IC AY I9N  T£9« 

08  1J0  J«l,Nl 
01  130  I»l, -1 

UC  93(  1  ,J )*83(I, jj*9,5* (>z U , J)*8v2( J)*v2 ( t , J*1 )*«vi 1 J*l ) J 
ENO 


0072000 
0C7J0C0 
007*000 
00 '5000 
007*000 
0077000 
00  78000 
0079000 
00*0000 
OOOIOOO 
00*2000 
00«30C0 
00**000 
00*5000 
00**000 
00*7000 
00*8000 
00*9000 
0090000 
00*1000 
00*2000 
00*3000 
oon*ooo 
O0*5«00 
009*000 
0097000 
00*8000 
00*9000 
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uioiooo 
o  i  r  2  c  o  c 
010  3OQ0 
oiomoqo 

01V500C 
OlOtCOO 
0107*00 
0108000 
0109000 
0110000 
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0112000 
0113000 
011*000 
0115000 
01 1*000 
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0 140000 
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•  **  “£**#£*  *uTOUT 

1. 

2* 

3* 

*• 

5. 

*• 

7* 

»* 

4a 

10a  C 
H*  C 
U*  C 
13* 

I«* 

IS* 

10* 

IT* 
l»* 

1*« 

JO*  1 

2»* 

22* 

23* 

jo* 

JS*  to 

JO* 

JT* 

JO* 

JO* 

30*  Jo 

Jl* 

32* 

33* 

3«. 

35*  3c 

30* 

JT* 

JO* 

jo. 

afl*  ac 

ala 

42* 

a  3* 

aa* 

«?*  5a 

ao* 

47* 

a6* 

«9« 


30«»0UTlNP  T 0 L T 

»*s»itTEa  -*21,n»<i 

0*R4»*tTEa  ,►?.«. 2 , N 1  a*..l , N  jtta.j 

COx^/ONE/VRI  (a,Ml,  VT  t  t  ** .  ►•  3  >  ,  VJt  t»l,R),Bl  ("1 , 4  1)  ,  V*2  <*,M  )  , 

C2-**9N/TH«/RH",R^eR<».n.8vj(to,*i.oH»,e».!i4,eNoo,C8»i,s,»«(-i,zIcti*) 

C0».»9*t/FflB/PE|.T,*M»E,lTlaiE#jSTE*,:8*<e,IT*PE,TBV 

T»ls  SueROUTIRE  RRJNT  BUT  7IEUDS  08»  *  OUlCr  LOOK 

riHtusi**.  iPu^t*,*) 

7*0  0*4*41  <////,'  44ei*l.  VEUCITY  CC*/3 }  *T  T*i,Ia,<  h') 

144CEMJ41  VELOCITY  CC'/S)  * T  Ta  1 , 1. ,  i  hi) 

VE4T3C4L  VELOCITY  fC*/S)  4T  T*1,!*,!  h  • ) 

EL0Y4NCY  HELD  (*.001)  iT  Ta<,l*,i  Hi) 

PRESSURE  (*I0  0YNE/C“**2 )  4T  T1*E»M»,'  hi) 


7"5  rSRNiT {////, i 

7 1  0  P*fi**l (////, • 

715  00p«»t (////, l 

725  0"R**1  . -s  ....  ... 

72c  0*R*4t(|h|#////#^ox# '♦♦♦♦♦♦♦*♦♦ 


0*.J*'  D4Y  ISTEO  4', 

u*YR«Ti*E/eaano,«t ,oce i 

PRINT  720,ITIME,D47,ISTEP 
0"  10  JR1,N1 
DO  10  1*1, R 
IOu**(l,J)»vRJ(I,J) 

PRINT  700.IT2-E 

C*UU  '’*P(!tl|i*«Rlf22»o,Nl) 

C"  JO  JR  1  (I1 1 
0*  Jo  1*1, “ 
lOyRd,  J)4VTJ(I,J) 

PpIfJT  705,IT|»f 

C‘U  *,4PnO'"*,R).22*’,#M) 
C*  JO  J*1,K 
r"  Jo  1*1, m 
I On“  tl  ,.’)*v  7j(  I  ,  J ) 

PRINT  71U.ITI-E 
C»lL  *4P ( loo* ,*2 • 2 1 (* 1 (N ) 
40  J«l,(,l 
C*  «.J  I«l,»l 
ICu-Ct,J)aBj(i(j).,,tj 

»°int  7i*. it s-e 

CHI  **P(IPip~,R«,2a,Yi  ,N,) 
O’!  5u  Jtl.lll 
0‘  5*  1*1,111 
irL,N(lfJl»P<t,J).t,£.| 
PRINT  J  5  » 1 T I  *1£ 

C»Lt  ’■»P<I0|.»'.Rc,2J,*1,'.U 

RE  TURN 

E*C 


17,' 


OUTPUT  4T  TI*E 
') 


•'iIO.i 


0* 


0001000 
OOOJOOO 
OOOJOOO 
000*000 
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000*000 
0007000 
OOUOOOO 
00O400C 
0010000 
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ooi  jooc 

0014000 
OOISOOC 
001*000 
0017000 
0010000 
00 1 4000 

oojoooo 
004 1000 
0042000 
004 JOOC 
0  0440  0  0 
004500C 
004*000 
0047000 
0042000 
1)04*000 
OPJOOOO 

ouJiooo 

ooJiooe 

0043000 

OcJaOOC 

ooisooo 

003*000 

0037000 

oc Jeooo 
coioooc 
00*0000 
oo“icoo 
00*2000 
no* jooo 
oc*«occ 

00*5000 

oc**ooo 

110*7000 

10*«('00 

00**000 
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*•  -t**t*»  “  AP 


1  • 

SUa-tlTISC  “lP(»,fc, 

oooirce 

!• 

paaa-fc Tt«  »6Jl,s»«l 

ooo jogo 

S« 

OI-CNSP*.  S<*a“),Z(t.K) 

0003040 

u  ■ 

JMEStl.  a 

OOOaiCOO 

5« 

TO 

r'la«tt(l>iS,?*,2*|«) 

oouseoo 

64 

»« 

0006000 

?• 

0007000 

¥  • 

p*  1(1 

0008000 

ta 

tfl 

!«(t)«»U)at.E.c*r,l 

0006000 

10a 

ao  J ■  t  , aj ** 

001000c 

1  la 

20 

tZ(J>*a(J)*t.E»c*C.l 

0011000 

ia* 

P»t*T  70 

ootaooo 

ti* 

P«t"T  TO,(IS(t),l«l,‘P) 

0013000 

i«* 

P»t’4T  TO 

0016000 

1  5a 

00  Jn  JJal.'m 

0015000 

16a 

JSNN* i«j j 

0016000 

l?» 

30 

ppjnt 

0017000 

1  6* 

«ETUa». 

OQISOOO 

|  6  • 

£*•0 

0Q190Q0 

C*"C  lHAGt  '1LE 


*»E  ®S  1  ***■  CS.J« 
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•  M  PRESS 


si ®«otiiNf-  P8£ji  oci'ioec 

retiecc 

f*-is  SL'S®»tT!KE  SETS  Ue  *0*CInG  FUNCTIONS  an’'  F0'lNf;ABY  00U  SC  V  C 
CONDITIONS  F»s  Txf  POfSSLHE  DIAGNOSTIC  EOLATIONS  EPS  SCOCUiCOO 

OpoSOQC 

pabaultEn  “*2t,N*;t  CCBpPOC 

parameter  Ml8“.|  ,  »2**-*2*M*'-«1  ,‘;2«n*2  O0O7000 

PARAMlTffc  N8l**2*NBL*1«NHL«*1  OOOFCOO 

REAL**  RC*P*RInv,SINV1,mTIlOA,C,U-hyI  00U4COC 

taxMON/Ol.F./VRl  (A  ,M  ),  v  T1  vz  1  (*1  »N),B1  (HI  ,M  >,VK?(«,*t ),  o«iocee 

!  vrj(-,KlJ»Vl2t«1,A>),B2C“l,Nt)f»8J(«#8J)«VTJ(M<AJ)f  OOUOOO 

2  V23(*M#K),BJ(**l,i.l),PtMl,N1)  0012000 

CAM»Qs/T.8sRt(««),«2i>M)(C*l("l),Di>2(“)(Zt(N!,Z2(M)«CZlfM)>DZ2(N)  OOJ  JOOO 
c,w»a*>/tH8/8R',FfeB  [M  J,BV2(n),alPhA,Bn0a,BNDR,CB«I*C,H.  (-)  ,Z*(N)  00  1«OOC 

0015000 

«■•«*.*  IS  Tm£  SIZE  "F  «  AND  F  00>*000 

►  P«®  l  «2  ,*iPiM*M(»L*«  1  001700C 

00X000 

PARAht  TER  »*!*'♦  |  ,AP*s«  I  00X000 

PARA-tTEB  hPiihP.J  ,  >>Pt«*,P-i,  no  1»».P.  1  ,KP2lNP«2  00*0000 

CPH-0N/Ev8/BJ>4v<Pf2,»,F2(NeL*)»«I'<V|  {•*F2,“F2#MBt*t  J,»C9B(*F»  J),  00*1000 


PABAMt  TER  »®BM» 1  ,ABlN* I  0014000 

PAAAMfcTEB  fPlBMF.)  ,»>P2»“P«f  ,NPl»l»».1  ,*.P2lNP-2  00*0000 

Cpp-0N/Ev»/Bjuv<pf 2,pF2,s8L«)/«Invi (**2 »“P2#NBl* 1 ) *®COS (*P , J) ,  00*1000 

1  RTILI  A(PP2),F(*.*,A,Pj,NBS122(>.BLA5#IS(t.PL«),SL,'FfSBL|0,  00*20  00 

2  tE(NfcL»  )»f  1 ICHP )  ,F)N  f«P  ),  F2 |(nP),F2H  (».P),A»CpF),  AT  <«P),  00*1000 

1  PB(“F,nP).CX(“P)»C*{N*1  OC*«000 

DIhENSIOn  OKAi'ATl  (»*2»'F2>,«<“F,t,p>  00*5000 

EQUIVALENCE  (CUHHY1  ,81 A  V  ( 1  ,  1  ,NRL*  )  J  00**000 

pabAmETEB  HPMPaPpAPP  00*7000 

DATA  i)U*‘"Vl/A>P«8*c(/(>/“8MP*0</  00*6000 

0 A T A  a*C*LL/0/  00*4000 

HB5IZ2  BtP«CSENTS  NUNBE*  OF  jrTebjpb  g*IC  FOISTS  I*  EACH  BLOC*  1“  X»DIOOiOOOO 
S2  BEPBESthTS  NUhR£a  OF  {KTERTO*  GRID  POINTS  IN  Y«DIRECtI8N  OoJlOOO 

NHL*  RE*RtSENTS  nu-nea  Of ' BLOCKS  IN  INDIRECTION  OOlSOOO 

THE  VARIABLES  At  1  ,A]N(a21  «A2H  TAaES  The  value  0  FOR  0IRICHL£T  B.C.  OOllOoO 

AND  1  fOH  nEUhAnh  t.c,  AT  THEIR  RESPECTIVE  BOUNDARIES  At,  CO*RES*BNOOO )«OO0 
Jal  At*.  TO  JBN  A2t  T"  tit  A2H  re  I»P  00)5000 

BOUNDARY  CONDITIONS  ARE  OOAPOOO 

*(I.U*U»AU)»X(I,I)*AII»(*(I,2)»*1|U))  00)7000 

*(J»NP)*O.AtN).I{I,NP)YA1N«{X(l,NP»OYFlN(I))  OOlSOOC 

l(l«J)P(l>A2l  l.*(l,J)«A2l  a(*(J,J).F2UJJ  J  00)4000 

*(HF,))#(l«A2H,.1,p#,J)«42H.(*{HP,t<J)pF2H(J))  00*0000 

nCALLPnCalL‘1  00*1000 

00*2000 

define  tnE  fopcjng  function  of  the  elliptic  equation  oo*jooo 

00*4000 

or  10  j*i,nP2  oo«5ooo 

c"  io  i*i,-*2  oo«*ooe 

10  F()*1#J41  latex CI •  1 1*0*211- *1  )«VR1(I*I ,  J)*VZ1C,  J*l  l/CZHj)  00*7000 

1  .A»U*t). 0*2(1). VRT(I,J).VZJCtfJ)/OZl(J))»RH*  00*6000 

00*4000 

set  lp  an  initial  guess  oosoooo 

0051000 
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M  605  JoJ, '<0  1 
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